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1 


Abstract 


In today’s society, there is a wide demand for high-precision and high-stability time service in the fields of electric 
power, communication, transportation and finance. At present, the time standard in various countries is mainly 
based on atomic clocks, but the frequency drift of atomic clocks will affect the long-term stability performance. 
Compared with atomic clocks, millisecond pulsars have better long-term stability and can complement with the 
excellent short-term stability of atomic clocks. In order to improve the long-term stability of the atomic timescale, 
and then improve the timing accuracy, this paper proposes an algorithm for steering the atomic clock ensemble 
(ACE) by ensemble pulsar time (EPT) based on digital phase locked loop (DPLL). First, the ACE and EPT are 
generated by the ALGOS algorithm, then the ACE is steered by EPT based on DPLL to calibrate the long-term 
frequency drift of the atomic clock, so that the generated steered atomic time follows both the short-term stability 
characteristics of ACE and the long-term stability characteristics of EPT, and finally, the steered atomic time is 
used to calibrate the local cesium clock. The experimental results show that the long-term stability of atomic time 
after steering is improved by 2 orders of magnitude compared with that before steering, and the daily drift of a local 
cesium clock after calibration is less than 9.47 ns in 3 yr, 3 orders of magnitude higher than that before calibration 
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1. Introduction 


High-precision timing is realized on the basis of establishing 
a high-accuracy and high-stability timescale, which is widely 
used in power, communication, transportation, finance and 
other fields. High-precision time reference is mainly generated 
by the joint use of atomic clocks, including common atomic 
clocks such as a rubidium clock, hydrogen clock or cesium 
clock, which have good short-term stability and small timing 
error in a short time of self-timekeeping. However, due to the 
long-term frequency drift of atomic clocks, in the scenario of 
long-term self-timekeeping, atomic clocks and timekeeping 
systems based on atomic clocks will produce more and more 
serious frequency drift and phase deviation, resulting in 
increasingly large 1 pulse per second (1PPS) timing errors. 
Coordinated Universal Time (UTC) is a widely used atomic 
time (AT), and the local timekeeping laboratory regularly sends 
the atomic clock comparison data of UTC Oh every day to the 
Bureau International des Poids et Mesures (BIPM). BIPM 
obtains the time difference between the local AT data and the 
international comparison center station through the remote time 
comparison link, and obtains the weighted average free AT. 
International Atomic Time (TAD) is calibrated with a fountain 
clock, and UTC is obtained by a leap second (Zhang et al. 
2020). 


Pulsars are high-speed rotating neutron stars with excellent 
long-term rotational stability. Pulsar signals are used to steer 
atomic clocks and generate timescales that take both long-term 
and short-term stability into account. The high-precision real- 
time observation of pulsar times of arrival (TOAs) by advanced 
telescopes such as FAST can avoid the lag problem of precise 
timekeeping sources such as UTC and Terrestrial Time at BIPM 
(TT(BIPM)), and can be further applied as timing signal sources 
in deep space navigation, power grid, intelligent transportation 
and other fields to achieve high long-term stability and high- 
precision timing (Luo et al. 2021; Qiu et al. 2022; Yang et al. 
2022). In addition, pulsar pulse signals are much more complex 
than global navigation satellite system (GNSS) timing systems 
that can be easily identified by machine learning methods (Tariq 
et al. 2022), and pulsars can serve as auxiliary signal sources if 
the existing atomic clock timing system is deliberately destroyed 
or tricked, or loses its lock in areas with weak signals (Qiu et al. 
2022). The generation of pulsar time (PT) is based on the pulsar 
timing observation data and ephemeris data, and the difference 
between the predicted pulse TOA and the observed value at the 
solar system’s centroid is obtained, that is, the timing residual. 
Here, the predicted value is the PT generated by the pulse phase 
model, and the observed value is based on local AT, so if all 
error terms are corrected, the timing residual is equal to the 
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difference between PT and AT (Zhou et al. 2021a, 2021b). By 
tracing the AT recorded by the atomic clock of the station to 
TAI, the clock error between PT and TAI is obtained. By 
combining the timing residuals of multiple pulsars with some 
algorithm, the ensemble pulsar time (EPT) can be established. 

Based on the characteristics of PT and AT, and the 
complementary characteristics of stability, the combination of 
PT and an atomic clock can improve the accuracy and stability 
of time reference. In recent years, in the main research on 
utilizing pulsars and atomic clocks for joint timekeeping or 
timing, Liu et al. (2023) proposed using Vondrak—Cepek 
filtering to combine PT and AT to generate a timescale with 
good short-term stability and long-term stability, but the 
timescale generated by the algorithm lags behind and the 
real-time performance is insufficient. Li et al. (2023) relied on 
the dual-steering algorithm to steer an atomic clock by pulsars, 
in which the monthly steer is closed-loop, which ensures the 
real-time performance of the time reference after steering 
(Dong 2020). However, the pulsar timing residual in this paper 
is a simulation of future data, and the longer the prediction 
time, the less credible the sequence is. 

In this paper, an algorithm based on the digital phase locked 
loop (DPLL) for steering the AT is proposed. The atomic 
clock signal is set as the local oscillator of DPLL, and the 
pulsar signal is set as the steering frequency source of DPLL. 
Most of the low-frequency noise of the atomic clock and most 
of the high-frequency noise of the pulsar are filtered by the 
loop filter, and a timescale combining the advantages of the 
two is generated, so as to generate a high-precision timing 
signal and provide 1PPS for timing sources such as power grid 
terminals. 

The parts of this paper are arranged as follows: Section 2 
introduces the principle of timekeeping algorithms, including 
the EPT algorithm and atomic clock ensemble (ACE) 
algorithm. Section 3 introduces the principle of pulsar steering 
algorithm based on DPLL. In Section 4, an experiment is 
carried out and the results are analyzed. The clock errors of 
EPT(1PPS), ACE(1PPS) and (ACE+EPT)(1PPS) based on 
DPLL are calculated, and a local cesium clock is calibrated 
with (ACE+EPT)(1PPS) to achieve high-precision daily drift, 
and the error between the experiment and the actual timing 
situation is analyzed. Section 5 summarizes this article. 


2. Principle of Timekeeping Algorithm 


This section introduces the basic principles of a timekeeping 
algorithm, including the EPT algorithm and ACE algorithm. 
Among those, the ALGOS algorithm is the basis of clock error 
weighting for pulsars and atomic clocks, and the corresponding 
1PPS clock error signals are generated by the EPT algorithm 
and ACE algorithm respectively, which are used as the external 
frequency source and local frequency source of DPLL. 
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The calculations of EPT and ACE need several pulsar clocks 
and atomic clocks respectively. When using these two clock 
groups to calculate the average timescale, methods such as 
ALGOS, Kalman filtering, and wavelet analysis can be applied. 
The ALGOS algorithm, as a typical weighted average 
algorithm, can improve the stability of EPT and ACE, so it is 
employed to calculate the comprehensive timescale. In this 
algorithm, the Allan variance used to calculate atomic clocks 
and the statistical method of o(7) variance (also known as 
Hadamard variance) used to calculate pulsar clocks are 
expressed as Equations (1) and (2) respectively (Allan 1966; 
Demetrios et al. 1997) 


2. — . = . +\72 
o\(T) = IN DT? pr 2 [x(i + 2) — 2x(i + 1)+ xF, 
(1) 
7 
aT) = (c). 2) 
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Among them, the residual and error are divided into 
subsequences with equal intervals of 7, and each subsequence 
is fitted by a cubic polynomial, where c3 is the coefficient of the 
highest term, and () is expressed in all subsequences, and the 
average value is taken as the weight inversely proportional to 
the square of c3 error. 

Pulsars have sizable, a priori unknown frequency drifts, so 
Allan variance is not particularly suitable for pulsar data 
(Demetrios et al. 1997). By defining one more difference from 
the Allan variance, o,(7T) variance ignores the fixed frequency 
drift in the pulsar timing residuals. 


2.1. Ensemble Pulsar Time Algorithm Flow 


In this section, first, the timing residual of the selected pulsar 
is used to generate EPT, which is traced back to TAI through 
TT(BIPM19) and converted to UTC, and then interpolated into 
EPT(1PPS) clock error signal as the steering frequency source 
of DPLL. The detailed flow chart is shown in Figure 1. 

In order to improve the accuracy of the final EPT(1PPS) 
timekeeping, according to the indexes of “more observation 
points” and “longer observation time” proposed by Shaifullah 
et al. (2022), M better pulsars are selected as frequency 
sources. Due to the strict linear relationship between PT 
and timing residual, we can use the Tempo2 software 
(Edwards et al. 2006; Hobbs et al. 2006, 2009) to fit the 
observation data and calculate the timing residual traced back 
to TT(BIPM19). 

Due to the large timing noise of a single pulsar, the timing 
stability can be improved by weighting. Therefore, it is 
necessary to set the reference pulsar to Modified Julian Date 
(MJD), and align the timing residuals and errors of other 
pulsars to the MJD reference by down-sampling. In order to 
calculate the clock error of EPT(1PPS) with EPT, it is 
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| Extract the original timing residuals of 9 pulsars. l 


Differentiate the EPT(1PPS) clock difference. | 


Figure 1. Calculation flow chart of EPT(1PPS). 


necessary to linearly interpolate the timing residuals and errors 
with the interval of 1PPS by MJD, so that the interpolation of 
subsequent timing residuals can meet the condition of equal 
interval and reduce the errors. 

Because the purpose of this paper is timing, it is necessary to 
align the past 1PPS signals of EPT and ACE and input them 
into the DPLL steering algorithm to fix their time series length. 
When adding the latest value, the earliest value is discarded. In 
order to calculate the long-term stability, the o(7) method 
should be used to evaluate the stability before truncation for 
weighting, so as to eliminate the frequency drift of the pulsar. 

Since the AT has been traced back to UTC, and the pulsar 
signal is traced back to TT(BIPM19), it is only necessary to 
trace the PT back to UTC, so that the time reference of the 
steering operation can be unified. The interval of traceability 
data between TT(BIPM19) and UTC provided by BIPM is 
10 days (Rodin 2008), with larger intervals, so it is necessary to 
linearly interpolate the timing residual points to calibrate the 
deviation from UTC. Then, we intercept the data according to 
the preset time series length. 

To facilitate interpolation, create an MJD index table such 
that the horizontal axis is a 1PPS time series, the minimum 
value of the vertical axis is O and the maximum value is the 
number N of timing residuals within the set MJD range. 
Distribute the abscissa of the 1PPS signal (1, 2, ...) evenly to 
the gaps in timing residuals. 

Let the abscissa of an interpolation point be x and the 
interpolation be y. If the abscissa of the original residual is 
replaced by x1, %2,...,.xy and the ordinate corresponds to yı, 
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| Pre-process: Interpolate to complete missing signals | 
of atomic clock measured data. | 


L Weighting 1PPS of 3 atomic clocks to generate ACE(1PPS) | 
clock error according to ALGOS algorithm. 


Figure 2. Calculation flow chart of ACE(1PPS). 


Y2,...,Yx respectively, the interpolation formula can be obtained 


yp I= — y), x<% 
y= Jp t ie OT Y) MH [Ê x S xy. (3) 
Yy + as ae = Yy_p. xX > xn. 


After the PT (1PPS) clock error of each pulsar is obtained by 
interpolation, the normalized weights w; are calculated based 
on the ALGOS algorithm, so that the clock error of M pulsars 
can be weighted to calculate the EPT(1PPS) clock error, and 
weaken the instability brought by a single pulsar (Petit et al. 
1993; Panfilo & Arias 2019). 

Finally, the EPT(1PPS) clock error needs to be differentiated 
to obtain the frequency difference, so as to complete the 
steering of the ACE(IPPS) phase in the subsequent DPLL 
steering algorithm. Because the pulsar signal is used as the 
steering frequency source, the steering physical quantity that 
needs to be provided is the frequency, not the phase. Since the 
EPT(1PPS) resulting from this step is a differential frequency 
signal, o(7) is no longer required for subsequent evaluation of 
its stability, but only Allan deviation. While the generated 
(ACE+EPT)(1PPS) has an even smaller frequency drift, Allan 
deviation can also be used to calculate stability. 

After the above steps, the real-time availability of the EPT 
signal is realized, and at the same time most of the original 
timing residuals of each pulsar are preserved, which are used as 
the steering frequency source of the subsequent DPLL algorithm. 


2.2. Atomic Clock Ensemble Algorithm Flow 


In this section, first, the measured clock error of the 
pretreated atomic clock is fitted and interpolated by a quadratic 
model, and the noise is added at the interpolation point and 
weighted based on Allan deviation, which is used as the local 
oscillation signal of DPLL. The detailed flow chart is shown in 
Figure 2. 

First, the actual clock errors of atomic clocks should be 
analyzed in detail, the clock error data with obvious abnormal 
values should be eliminated and linear interpolation should be 
performed for very few clock error missing points. If the 
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original data are normal, skip this step. This step can ensure the 
normal progress of subsequent data processing, and at the same 
time, due to the small number of abnormal data points, the 
impact on the original data can be ignored. 

If the time interval of the original data is too long, it needs to 
be interpolated into 1PPS clock error. Generally, atomic clocks 
have clock drift that deteriorates with time, so the quadratic 
clock error model (Galleani et al. 2003) is often used to model 
them, and the approximation of clock error, clock speed and 
clock drift parameters is completed. The expression is as follows 


1 
x(t) = xo + Yot + al + e,(t). (4) 


Here xo, yo and d are initial time error, frequency error and 
linear frequency drift, respectively, which are deterministic 
components, and the rest are noise. 

In order to steer the atomic clock by pulsars, it is necessary 
to align the MJD of the selected atomic clock with the MJD of 
EPT(1PPS) clock error by interception, and interpolate the 
clock error of the atomic clock into 1PPS clock error with the 
fitted parameters. Similar to the data processing method of a 
pulsar, Allan deviation (Allan 1966) can be calculated 
according to historical observation data before intercepting 
the atomic clock error, so as to calculate the more accurate 
weight of each atomic clock. If the 1PPS signal can be obtained 
directly in the original observation, this step can be skipped. 

Based on the ALGOS algorithm, the weight of each clock 
can be obtained by Allan deviation after generating the 
AT(IPPS) clock error, and then the ACE(1PPS) clock error 
can be generated by weighting. After the above steps, the 1PPS 
matching between the ACE signal and the EPT signal is 
realized, and it is used as the local oscillator frequency source 
of the subsequent DPLL algorithm. 


3. Principle of Pulsar-steering Algorithm Based 
on DPLL 


This section introduces the algorithm principle from time- 
keeping to timing. The timekeeping algorithm completes the 
steer of EPT(1PPS) to ACE(1PPS) based on DPLL, in which 
the steering frequency source and local oscillator frequency 
source of DPLL are static. The timing algorithm is based on 
(ACE+EPT)(1PPS) of steering output, which periodically 
calibrates the clock error of a local cesium clock and realizes 
the dynamic output of actual physical signals. 


3.1. Principle of Loop Establishment and Steering 
of DPLL 


DPLL is a negative feedback system, which can steer the 
phase of a local oscillator through the phase of the input signal 
source. Its basic structure is illustrated in Figure 3 (Prasad & 
Sharma 2012). 

Among that, PD, LPF and VCO are phase detectors, loop 
filters and voltage-controlled oscillators, respectively. The 
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(ACE+EPT) 
(1PPS) 


EPT (1PPS) 
(1PPS) 


Figure 3. DPLL schematic diagram. 


product of the gains of the three is recorded as G(z), which 
represents the transfer function of the DPLL open-loop system in 
the Z domain; EPT(1PPS) stands for the reference input of 
DPLL, also known as the steering frequency source; ACE 
(1PPS) represents the local oscillator frequency source of DPLL, 
input to the VCO; (ACE+EPT)(1PPS) is the output after 
steering, fed back to the input terminal through the frequency 
divider. By analyzing the J/O characteristics of the DPLL 
system in Z domain, the output of DPLL can be expressed as 


G(z) 
ACE + EPT = EPT } ACE 
( + ) (z) iGO (z) 1+ 6® (z) 
= H(z)EPT(z) + He(z) ACE(z). 
(5) 


Here H(z) is the transfer function of the closed-loop system, 
and it is a low-pass filter, which can filter out the high- 
frequency noise in the clock error of the steering frequency 
source EPT(1PPS). H(z) is the closed-loop error transfer 
function of DPLL, and it is a high-pass filter, which can filter 
out the low-frequency noise in the clock error of local oscillator 
frequency source ACE(1PPS). Therefore, the steered (ACE 
+EPT)(1PPS) clock error has small low frequency noise and 
high frequency noise at the same time. 

Considering that the third-order phase-locked loop can 
stably track the frequency ramp-up signal (Lu 2016), and its 
performance is better than that of the second-order phase- 
locked loop, which can ensure the steering stability of ACE 
(IPPS) with EPT(1PPS) clock error in this paper, the system 
transfer function of DPLL is selected from Feng et al. (2021) 


A(z)= 


KıKo (2# T z! ) 
+ 1 


1— z! 271 j71-—z 


-1 
, T 1l+z i (6) 
In -T1-—z! 


Here T is the sampling period, K, is the gain of the phase 
detector, Ko is the gain of the voltage steered oscillator and the 
intersection bandwidth of the system transfer function is 
_ Ka kore 


277} 


K (7) 
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The optimal loop parameters (Shen 2018; Feng et al. 2021) are 
obtained when and only when the single sideband phase noise 
intersection of ACE(1PPS) and EPT(1PPS) is equal to K. At 
this time, (ACE+EPT)(1PPS) simultaneously retains the near- 
end phase noise of EPT(1PPS) and the far-end phase noise of 
ACE(1PPS), and the steering effect is the best. 

In the actual timing scenario, according to the dynamic 
change of input sequence, the value of K is automatically 
adjusted according to the calculation formula of loop parameter, 
and the loop can remain stable. When the output signal is in the 
first steering cycle, the 1PPS clock error of the steered source 
ACE is output. In the following steering cycle, the frequency of 
ACE(IPPS) is feedback steered by calculating the difference 
between the frequency of EPT(1PPS) and that of ACE(1PPS), 
thus avoiding the uncontrollable phase drift of ACE(1PPS). At 
the same time, using the excellent short-term stability of ACE 
(PPS) clock error, a lot of high-frequency noise of EPT(1PPS) 
clock error is eliminated, and the second-order difference of the 
output signal is similar to that of ACE(1PPS) clock error. 


3.2. Daily Drift Calibration of Local Cesium Clock Error 


In the last section, the (ACE+EPT)(1PPS) clock error signal 
in the past year was obtained by DPLL. However, only the 
1PPS clock error in the last second of this output signal can be 
used for timing, because only the 1PPS clock error in the last 
second of the two input signals is real-time. Although the clock 
error of an atomic clock is better predicted according to the 
model, pulsars are affected by various complex factors, and the 
uncertainty of clock error prediction is too strong (Demetrios 
et al. 1997), so it is difficult to accurately generate future 1PPS 
clock error signals. Therefore, the clock error of (ACE+EPT) 
(1PPS) in the last second is selected, and then the one with the 
smallest daily drift among L atomic clocks is selected, and the 
(ACE+EPT)(1PPS) clock error is used as its calibration source 
to obtain CSgteereq(1 PPS). 

When the DPLL runs once a day, first, the ACE(1PPS) clock 
error and EPT(1PPS) clock error of that day should be 
introduced in the same way, and the clock error of the earliest 
day should be deleted. Then the initial phase of the current one- 
year sequence is calibrated to UTC released one year lagged. 
Finally, the 1PPS clock error in the last second is calculated by 
the same method. At this time, the local master clock is 
calibrated by the method of phase microjump (Fang et al. 2022) 
to avoid the loss of lock and the deterioration of stability 
caused by phase step. Therefore, the final daily clock error is 
the sum of the daily phase offset and calibration error of the 
local atomic clock. 


3.3. Algorithm Flow of Pulsar Steering Atomic Clock 
Based on DPLL 


The flow chart of clock error steering algorithm based on 
DPLL is depicted in Figure 4. 
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Figure 4. Flow chart of ACE(1PPS) calculation. 


These include EPT algorithm, ACE algorithm and clock 
error steering algorithm based on DPLL. The final generated 
CSsteered(1PPS) is the output signal of this algorithm, which is 
used for high-precision timing. The clock error steering 
algorithm based on DPLL is as follows. 

Step 1. Calculation of ACE(1PPS) and EPT(1PPS); 

Step 2. Based on the steering of DPLL clock error, the 
timekeeping signal of (ACE+EPT)(1PPS) is generated; 

Step 3. Calibrate the local cesium clock periodically with 
(ACE+EPT)(1PPS) to complete high-precision time service. 

To sum up, the algorithm in this paper can be summarized 
as: EPT timekeeping algorithm represented by red dotted line 
box, ACE timekeeping algorithm represented by green dotted 
line box and DPLL-based clock error control algorithm 
represented by blue dotted line box. Although the algorithm 
is relatively complicated, it can realize automatic processing in 
the actual time service process, which can generate real-time 
time service scale with better accuracy and long-term stability 
than atomic clocks. 


4. Experiments and Results Analysis 


This section first introduces the source and selection of the 
data used, and then gives the calculation and analysis of EPT 
(1PPS) frequency error and ACE(1PPS) clock error respec- 
tively. Finally, the EPT(1PPS) frequency error is used to steer 
the ACE(IPPS) clock error through DPLL, and the local 
cesium clock is calibrated daily with the generated signal, and 
the robustness of loop parameter is analyzed. In each data 
processing link, the error compared with the actual timing 
system is analyzed. 
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Table 1 Table 2 

Key Parameters of Nine Pulsars (Arzoumanian et al. 2015) Key Parameters of Three Atomic Clocks 
PSR Name Span/yr TOAs MJD Range (days) Sampling Data 
30437-4715 13.50 5302 51770-56705 tome Sampling Period MID Missing 
J0613-0200 16.05 9322 50932-56796 Clock Name Points (days) Range (days) Ratio 
J1012+5307 16.83 13056 50647-56795 Cs1350104 360 5 53739-55534 0% 
J1643-1224 20.14 8136 49422-56778 Cs1351463 360 5 53739-55584 2.78% 
J1713+0747 22.46 17486 48850-57052 H1400711 360 5 50647-55564 1.67% 
J1744-1134 19.88 9834 49729-56991 
J1909-3744 10.82 11483 53041-56992 
J1918-0642 12.80 9941 52095-56769 
J2145-0750 19.83 8455 49518-56761 


4.1. Experimental Data 


Pulsar timing data are selected from the pulsar data set 
(Perera et al. 2019; Antoniadis et al. 2022) published by IPTA 
in 2019. According to Shaifullah et al. (2022), JO437-4715, 
J0613-0200, J1012+5307, J1643-1224, J1713+0747, J1744- 
1134, J1909-3744, J1918-0642 and J2145-0750 are selected 
from 65 pulsars. The key parameters of the selected pulsar are 
shown in Table 1. 

Among them, the observation range and TOA data points of 
the above pulsars are better in the pulsar data set, and the MJD 
range can cover the selected atomic clock data. 

AT data were obtained from the BIPM official website 
(https: / /webtai.bipm.org /ftp/pub/tai/data/), and three atomic 
clocks (Cs1350104, Cs1351463 and H1400711) from the 
United States Naval Observatory were selected (Liu et al. 
2023). The MJD range was 53739-55584, and the sampling 
period of clock error was 5 days. Among them, the data of 
cesium clock Cs1350104 are normal, while the other two 
clocks are missing less than 3% of the data. It is necessary to 
linearly interpolate these data, and then add noise with the 
average value of the second difference of the clock error as the 
average value, and then calculate the Allan deviation of each 
atomic clock error. The key parameters of the selected atomic 
clock are shown in Table 2. 

Because the DPLL steering algorithm requires two input 
frequency sources to have the same start and end time and the 
same data length, it is necessary to intercept the pulsar and atomic 
clock data with the same MJD respectively when verifying 
the feasibility of the algorithm. The ranges of three groups of 
MJD used in this paper are 53739-54104, 54104-54469 and 
54469-54834, respectively, and the length is one year. 


4.2. Experiment and Result Analysis 
4.2.1. Performance Analysis of EPT(1PPS) Clock Error 


This section analyzes the performance of EPT. 

Considering that the release time of UTC is one year behind, 
the timescale of 1PPS clock error is set to one year, so that the 
initial phase of DPLL local oscillator is aligned with the latest 


UTC. Since the AT data start from MJD 53739, in order to 
ensure a real-time and causal system, the residual of each pulsar 
after MJD 54104 is removed. J1909-3744 with the least 
number of residual points is selected as the benchmark of MJD 
to ensure the authenticity of timing residuals. Align the MJD of 
the other eight pulsars to the benchmark, then downsample and 
linearly interpolate the timing residuals and errors, then 
calculate the o,(r) variance of each pulsar. To maximize the 
use of the properties of the pulsar’s long-term stability, the EPT 
is weighted by the inverse of the end value of o(7). The o,(7) 
variance of nine pulsars and EPT are displayed in Figure 5(a). 

From the curve of EPT, it can be seen that the short-term 
stability and long-term stability of EPT are close to the pulsar 
with the best performance, especially on the one-year scale, and 
its o,(7) stability is better than all pulsars. Therefore, weighting 
the timing residual can smooth the uncorrelated timing noise of 
multiple pulsars and make EPT more stable than what is 
allowed by a single pulsar. At the same time, setting the length 
of pulsar timing residual to one year can ensure smaller timing 
error and better stability. 

IPTA Data Release 2 (DR2), the pulsar data set selected in 
this paper, has been traced back to TT(BIPM19). The 
calibration deviation of tracing is expressed as At(yus), then 
we have 

a = TT(BIPM19) — TAI — 32.184 s (8) 
UTC = TAI — 37s 


so that the PT can be traced back to UTC. Subsequently, the 
traced timing residual between MJD 53739 and MJD 54104 is 
intercepted. 

Based on MJD benchmark pulsar J1909-3744, the number of 
points for the timing residual is 1060. We substitute N = 1060 
into Equation (3), interpolate the timing residual to generate 
PT(IPPS) clock error and then add the evenly distributed 
random noise. We then substitute M=9 into the ALGOS 
algorithm to calculate the EPT(1PPS) clock error and calculate 
the frequency error through difference. The clock error of EPT 
(1PPS) is shown in Figure 5(b). From this figure, we can see 
that compared with the clock drift phenomenon of atomic 
clocks, the clock drift of pulsars is not obvious, and the UTC 
traceability error on a one-year scale is still in the microsecond 
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Figure 5. o(7) variance and EPT clock error of pulsars before and after weighting. (a) o,(7) variance of nine pulsars and EPT. (b) EPT clock error. 


order, indicating that the long-term clock drift of atomic clocks 
can be calibrated by pulsar signals. 

To sum up, the stability and clock error of EPT(1PPS) are 
better than those of a single pulsar. 

When actually calculating the EPT(1PPS) clock error, there 
are differences in data processing methods or systematic errors 
in the following steps: 


(1) When extracting the timing residuals of pulsars, it is 
necessary to use an astronomical telescope to collect the 
TOA and ephemeris data in real time, and generate the .tim 
file and .par file to calculate the timing residual (Edwards 
et al. 2006). Because the pulsar observation accuracy of 
the FAST facility has surpassed the data source IPTA 
(Yang et al. 2023) in this paper, the accuracy of timing 
residuals can be improved by establishing cooperation 
with similar observation stations. 

(2) When the reference star is selected by MJD and the other 
stars are down-sampled, the signal acquisition interval of 
the IPTA data set is uneven, the frequency band coverage 
is incomplete and the observation accuracy is limited 
(Perera et al. 2019), but in reality, the telescope 
continuously collects signals in real time, so the influence 
of this step on the timing residual will be reduced. 

(3) When tracing the residual to UTC through TT(BIPM19) 
and truncation, due to the inability to obtain TT(BIPM23) 
in real time, it can only be based on the calculation 
formula in TT(BIPM21) 


TT(BIPM21) = TAI + 32.184 s + 27667.5 ns 


—0.01(MJD—59579)ns = 0) 


to trace TT to TAI and further trace it to UTC. Since the 
difference between the predicted value and the actual 
value of TT (BIPM21) is between —6.2 ns and —0.2 ns, 


Table 3 
Deterministic Parameter Fitting of Three Atomic Clocks 
Atomic Clock Name Xo Yo d 
Cs1350104 4.1146 x 107 6.7652 x 10 6.2721 x 1071! 
Cs1351463 1.7747 x 107 4.4457 x 1078 3.7925 x 1071! 
H1400711 3.7494 x 104 9.9961 x 107 7.6749 x 107° 


the error will increase slightly in the actual timing system, 
but it is stable and controllable. 

(4) When linearly interpolating the residual and generating 
the PT(1PPS) clock error, the traceability error is less 
than 1%, corresponding to the original timing residual, so 
the interpolation effect brought by the sampling uni- 
formity of the timing residual will be more significant. 

(5) When generating the EPT(1PPS) clock error based on the 
residual stability weighting, the calculation error of the 
o(T) weight will be reduced because the TOA interval of 
the real-time observation of the pulsar array is more 
uniform, which makes the calculation more accurate. 


Generally speaking, the clock error calculation of EPT 
(1PPS) can be more accurate through the above steps. 


4.2.2. Performance Analysis of ACE(1PPS) Clock Error 


This section analyzes the performance of ACE. 

According to the quadratic clock error model of 
Equation (4), the clock errors are fitted, and the results are 
shown in Table 3, with complete clock interpolation based on 
the fitting parameters. 

Because this paper mainly uses the short-term stability of 
atomic clocks, the starting point of Allan deviation of each 
clock is substituted into the ALGOS algorithm, and the 
normalized weights of Cs1350104, Cs1351463, and H1400711 
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Figure 6. Clock errors between three atomic clocks and ACE(1PPS) and their Allan deviations. (a) Clock errors. (b) Allan deviations. 


are obtained as follows: 0.0509, 0.0356 and 0.9135 respec- 
tively. After generating the ACE(IPPS) clock error by 
weighting, we can plot it together with the clock errors of 
three atomic clocks as shown in Figure 6(a), and display the 
ACE(1PPS) clock error and the Allan deviation of three atomic 
clocks together as featured in Figure 6(b). 

In Figure 6(a), due to the high weight of H1400711, the 
clock error of ACE(1PPS) is close to it, and the relative UTC 
deviation of ACE is smaller than that of a cesium clock in the 
short term and hydrogen clock in the long term. The secondary 
clock error model is still dominant in the ACE, but the random 
noise term contained in it determines the magnitude of the 
Allan deviation diagram curve shown in Figure 6(b), and it can 
be concluded that the short-term stability of the ACE is better 
than that of the cesium clock. Moreover, the long-term stability 
is better than that of a hydrogen clock, and it is better than that 
of a single atomic clock as a whole. 

When actually calculating the ACE(1PPS) clock error, there 
are differences in data processing methods or systematic errors 
in the following steps: 


(1) In the selection of the measured clock error of atomic 
clocks, data from the National Time Service Center 
(NTSC) can be used instead (Zhang et al. 2020), and the 
time service stability and system robustness can be 
improved by increasing the number of atomic clocks. At 
the same time, NTSC has more public timing data, and 
even if the 1PPS signal cannot be read directly, the 1PPS 
clock error of each clock can be approximately solved by 
matrix operation. Therefore, the error in the actual timing 
system will be reduced and the system will be more 
robust. 


Table 4 
Power Law Spectrum Fitting Results of ACE(1PPS) and EPT(1PPS) 
Frequency 
Source hy ho hı h2 h4 
ACE(1PPS) 5x 10-7 0 0 5x10 3x 10° 
EPT(1PPS) 1x107 3x10” 0 0 0 


(2) When the quadratic model is used to fit the clock error 
and interpolation, if the 1PPS signal can be read directly, 
this step can be skipped, thus avoiding the simulation 
error in the statistical sense. Otherwise, the algorithm 
steps and errors are unchanged. 


Generally speaking, the calculation of clock error of ACE 
(IPPS) in practical engineering will be more accurate than the 
results in this paper. 


4.2.3. DPLL Steering Performance Analysis 


This section analyzes the performance of DPLL-based PT 
steering AT. First, the Allan deviations of ACE(1PPS) and EPT 
(PPS) are fitted separately based on the power-law spectrum 
parameters. The fitting parameters are shown in Table 4, and 
the fitting curves are displayed in Figure 7(a) and (b). Second, 
the single sideband phase noise of these two frequency sources 
is calculated based on the power law spectrum parameters, as 
shown in Figure 7(c). We calculate the intersection of phase 
noise to obtain loop parameter K. Finally, based on this 
parameter, the clock errors of ACE(1PPS), EPT(1PPS), and 
(ACE+EPT)(1PPS) before and after DPLL steering algorithm 
are drawn in Figure 7(d), where the EPT(1PPS) clock error 
represented by the green line is the pulsar clock error curve in 
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Figure 7. Loop parameter fitting and clock error before and after steering. (a) Allan deviation fitting result of ACE(1PPS). (b) Allan deviation fitting result of EPT 
(IPPS). (c) Single sideband (SSB) phase noise of ACE and EPT. (d) Clock errors before and after EPT steering ACE. 


Figure 5(b). The short-term stability of the blue line and red 
line is good, and their second-order difference is further 
compared and verified in Figure 8(a)-(c). 

From the fitting results of power law spectrum parameters 
shown in Table 4, it can be known that there is a small amount 
of frequency random running noise (h4) in the ACE(1PPS) 
clock error, which will greatly reduce the long-term stability of 
the atomic clock. However, EPT(1PPS) clock error only has 
phase white noise (2) and frequency white noise (họ), so its 
long-term stability is much better than ACE(1PPS). 

Since Figure 7(a) and (b) complete the accurate fitting of the 
frequency source, the single sideband phase noise calculated by 
the power-law spectrum parameters of ACE(1PPS) and EPT 
(1PPS) is accurate. By taking the intersection of the two curves 
in Figure 7(c), we can get 


K = 4.9091 x 10-°. (10) 


It can be seen from Figure 7(d) that the (ACE+EPT)(1PPS) 
clock error shown in the red line overcomes the clock drift 
defect of ACE, reduces the short-term clock error fluctuation of 
EPT and realizes the steering of ACE(1PPS) by EPT(1PPS). 
In Figure 8, Figure 8(a), representing the second-order 
difference of ACE(1PPS) clock error, is close to Figure 8(b) 
representing the second-order difference of (ACE+EPT) 


(IPPS) clock error, which demonstrates that the DPLL 
algorithm keeps the high-frequency part of the former well. 
Figure 8(d) is obtained by the difference between them. The 
clock error is reduced by about 2 orders of magnitude, and the 
curve trend is similar to the second-order difference of the 
frequency difference of EPT(1PPS), which shows that the 
algorithm keeps the low-frequency part of the clock error of 
EPT(1PPS) at the same time. 

We can enlarge Figure 8(b) as shown in Figure 9(a), then 
calculate the Allan deviation of ACE(1PPS) clock error, EPT 
(IPPS) frequency difference and (ACE+EPT)(1PPS) clock 
error, as depicted in Figure 9(b). 

In Figure 9(a), the clock error of (ACE+EPT)(1PPS) in one 
year is less than 3.60 ns. In Figure 9(b), the Allan deviation of 
(ACE+EPT)(1PPS) clock error first follows ACE(1PPS), then 
is biased by EPT(1PPS) and always follows its slope, and the 
stability after steering for one day is better than that of two 
signal sources. As the steering time continues to grow, (ACE 
+EPT)(1PPS) becomes more and more stable, and the one-year 
stability is 2 orders of magnitude better than that of ACE(1PPS) 
before steering. 

Although Equation (3) makes it difficult for the interpolated 
data of EPT to maintain the short-term stability of the original 
residual, resulting in the short-term Allan variance of EPT to 
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Figure 8. Second-order difference comparison before and after DPLL steering. (a) Second-order difference of ACE(1PPS) clock error. (b) Second-order difference of 
(ACE+EPT)(1PPS) clock error. (c) Second-order difference of EPT(1PPS) frequency difference. (d) (ACE+EPT)(1PPS) clock error minus the second-order 


difference of ACE. 


x10 


Clock error to UTC(s) 


EPT+ACE 


1.5 
Time 7(s) 


x10" 


Allan deviation a7) 


108 


101 


102 


=ð 

oO 
Es 
P 


ACE 
EPT 
EPT+ACE 


1016 


10° 


104 
Time 7(s) 


108 108 


Figure 9. Clock error and Allan deviation of MJD 53739-54104. (a) Clock error after EPT steering ACE. (b) Allan deviation before and after EPT steering ACE. 


10 


Research in Astronomy and Astrophysics, 24:035019 (12pp), 2024 March 


a 


x10? 
3 F T 


EPT+ACE 


Clock error to UTC(s) 


15 2 
Time 7(s) 


(e) 
D 
a 
a 


-9 | (c) 


Clock error to UTC(s) 


1.5 2 2.5 3 
Time 7(s) 


(=; 
>| 
a 
oe 


Allan deviation an) 


Zheng et al. 


ACE 
EPT 
EPT+ACE| | 


19°76 L ` L 
10° 10? 104 10° 108 
Time 7(s) 
(d) 
108 r 
ACE 
EPT 

T 4010 EPT+ACE | | 
5 
= 
S 
w -12 
S 10 
© 
pei 
E 
is} 
T 4074F 

40°16 f n f 

10° 10° 104 10° 108 
Time 7(s) 


Figure 10. Clock error and Allan deviation in different years. (a) The steered clock error of MJD 54104-54469. (b) The steered Allan deviation of MJD 54104-54469. 
(c) The steered clock error of MJD 54469-54834. (d) The steered Allan deviation of MJD 54469-54834. 


reflect the true value, the Allan variance of (ACE+EPT)(1PPS) 
follows ACE(IPPS) on a short-term scale and the clock 
difference of ACE(1PPS) over the steering period. The effect of 
imprecise and poor Allan variance of EPT can be ignored. 


4.2.4. Performance Analysis of Local Cesium Clock 
Calibration 


In this section, based on the (ACE+EPT)(1PPS) of the 
DPLL steering output, the performance analysis of the clock 
calibration of the local atomic clock is carried out, so that the 
steering output has a physical clock carrier and can be applied 
to the actual timing system. Since the total error of the timing 
system built in this paper is equal to the steering error plus the 
error caused by the local atomic clock self-timekeeping, by 
selecting an atomic clock with the smallest daily drift and 
performing a DPLL steering operation every single day, the 
1PPS timing signal with the best performance can be generated. 
For the steering error of 3.60 ns in Figure 9(a) (this value is the 
maximum value on the vertical axis), the cesium clock 
Cs1351463 with the slowest clock accumulation among the 
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three atomic clocks (it produces a maximum clock error of 4.73 
ns in one day) is selected, and the clock error and the control 
error are added together, we can see that the daily offset 
Of CSsteerea( 1 PPS) is less than 8.33ns in one year. 

Two additional independent experiments were performed 
below for added confidence. Considering that the timing residual 
of pulsars is the biggest variable in the steering experiment, the 
MJD of nine pulsars was directly replaced by 54104-54469 and 
54469-54834 in repeated experiments, and the loop parameters 
were unchanged. We summarize all the variables involved in 
atomic clocks and pulsars in different years into Table 5, and 
redraw the clock error and Allan deviation after EPT steering 
ACE respectively, as shown in Figure 10. 

As can be seen from Table 5, the weight of the atomic clock 
is slightly different in different 3 yr timespans, because the 
parameters of the atomic clock model are relatively fixed but 
the noise is random. The main pulsar and the number of 
residual points are almost unchanged, which demonstrate the 
system’s robustness brought by MJD homogenization. As can 
be seen from Figure 10(a) and (c), the maximum error of (ACE 
+EPT)(1PPS) is 4.74 ns, so the daily drift of Cssteerea( 1PPS) is 
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Table 5 
Comparison of Variables between Atomic Clock and Pulsar in Different Years 


MJD Range(days) Cs1350104’s Weight 


Cs1351463’s Weight 


H1400711’s Weight Main Pulsar Residual Points 


53739-54104 0.0508664660822 0.0356065262575 
54104-54469 0.0508664660801 0.0356065262561 
54469-54834 0.05086646608 10 0.0356065262567 


less than 9.47 ns in 3 yr, which is 3 orders of magnitude higher 
than that of microsecond-level daily drift caused by cesium 
clock independent timing before calibration. As can be 
ascertained from Figure 10(b) and (d), the Allan deviation of 
the output after replacing the input sequence still follows the 
short-term stability of ACE(1PPS) and the long-term stability 
of EPT(1PPS) steadily, which affirms that the steering based on 
DPLL is insensitive to the initial loop parameter and the system 
is robust. 

To sum up, the annual stability of AT based on the DPLL 
steering algorithm is improved by 2 orders of magnitude, and 
the daily drift accuracy of a local cesium clock calibrated by 
this signal is improved by 3 orders of magnitude. In the actual 
time service system, the clock calibration of the local cesium 
clock cannot jump, so the clock calibration is quadratic, 
compensated by clock speed and clock drift, and realizes the 
calibration completely after a period of time. According to the 
running time of the algorithm, the clock error can be calibrated 
by calibrating the clock reading at the previous moment, and 
also by using the method of phase microjump. In addition, the 
data length and calibration period can be slowly adjusted with 
the input data to maintain a high-precision time reference. 


5. Conclusion 


In order to improve the long-term stability and timing 
accuracy of AT, this paper proposes a pulsar-time-steering 
algorithm based on DPLL. Pulsar signal and atomic clock 
signal are used as the steering frequency source and local 
oscillator frequency source of DPLL, respectively. The finally 
generated steered frequency source not only avoids the 
frequency drift of the atomic clock signal, but also reduces 
the large frequency jitter of pulsar timing residuals in the short 
term. The experimental results show that, based on the 
automatic adjustment of loop parameters, the long-term 
stability of AT after steering follows the signal of EPT 
(1PPS), and the daily drift of (ACE+EPT)(1PPS) is less than 
4.74 ns in 3 yr, and the timing accuracy is improved by 3 orders 
of magnitude compared with that without calibration. After 
(ACE+EPT)(1PPS) is used to calibrate the local cesium clock, 
the daily drift of the clock is less than 9.47 ns. By applying the 
algorithm used in this paper to atomic clocks with smaller 
power-law spectrum noise, or introducing the low-frequency 
filtering method to generate EPT with higher long-term 
stability (Pizzocaro et al. 2020), higher timing accuracy can 
be obtained. 
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0.9135270076602 J1909-3744 1060 
0.9135270076639 J1909-3744 1060 
0.9135270076624 J1909-3744 1059 


The algorithm in this paper can be applied to power grid 
timing and other fields that need high precision and high 
stability. In the next step, we only need to replace the historical 
data used in this paper with the real-time atomic clock error 
(updated every five days, for example) and the real-time pulsar 
timing residual, so that we can complete the real-time calibration 
of the local atomic clock and improve the timing accuracy. 
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